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Since the kinetic and the potential energy term of the real time nonlinear Schrodinger equation can 
each be solved exactly, the entire equation can be solved to any order via splitting algorithms. We 
verified the fourth order convergence of some well known algorithms by solving the Gross-Pitaevskii 
equation numerically. All such splitting algorithms suffer from a latent numerical instability even 
when the total energy is very well conserved. A detail error analysis reveals that the noise, or ele- 
mentary excitations of the nonlinear Schrodinger, obeys the Bogoliubov spectrum and the instability 
is due to the exponential growth of high wave number noises caused by the splitting process. For a 
continuum wave function, this instability is unavoidable no matter how small the time step. For a 
discrete wave function, the instability can be avoided only for At k^ nax ^2TT, where k ma x = ir/Ax. 
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I. INTRODUCTION 



Taha and Ablowitai have shown for some time that the first order pseudo-spectral, split-operator method is a 
Qh very fast way of solving the nonlinear Schrodinger equation. Bandrauk and Shen- later applied higher-order splitting 
algorithms with negative coefficients to solve the same equation. They regarded the nonlinear potential as time- 
dependent. Since they can only estimate the intermediate-time nonlinear potential to second order, it is not proven 
that their higher-order algorithms actually converge at fourth or sixth-order. Recently Javanainen and Ruostekoski 3 
have shown by symbolic calculations that fourth-order algorithms are possible by use of the "latest" intermediate wave 
function in evaluating the nonlinear potential. Straucb^, by constructing a special operator that correctly propagates 
the nonlinear potential term, proved that this use of the "latest" intermediate wave function is valid. 

This work shows that: 1) Javanainen and Ruostekoski's finding is a direct consequence of Taha and Ablowitz' 
original work and a much simpler proof than that of Strauch is possible. 2) The time-dependent potential method of 
Bandrauk-Shen and the time-independent approach suggested by Javanainen and Ruostekoski both yielded identical 
O I second-order algorithms but different higher-order algorithms. 3) Verified numerically that algorithms derived by 
l— ~~ '■ the time- independent method do converge to fourth-order when solving the Gross-Pitaevskii equation. 4) All such 
splitting algorithms possess a latent numerical instability which causes the wave function to blow up despite excellent 
, total energy conservation. 5) The instability is shown to be due to the exponential growth of high wave number noises 
\Q ' intrinsic to the splitting process. For a continuum wave function, this instability is unavoidable no matter how small 
Q\ , is the time step. For a discrete wave function, this can only be avoided if At^2n/kf nax , which forces At to be very 
small if the discretization is very fine with a large k mas = n/ Ax. The next three sections summarize how higher order 
| algorithms can be systematically derived and Section V discusses the instability in detail. 

d 

fx ! II. SOLVING THE NONLINEAR SCHRODINGER EQUATION 

o 



Consider the nonlinear Schrodinger equation defined by 



*^ = (-2 V (2-1) 
The free particle propagation can be solved exactly in operator form 

<P(At) = e~ lAtf ?/>(0) (2.2) 

where the operator T = — ^V 2 . Since T is diagonal in k-space, (|2.2p is usually solved by Fast Fourier Transforms 
(FFT). Surprisingly, as shown by Taha and Ablowitz, the potential part of the equation 

4^ = <?MV- (2.3) 

can also be solved exactly 

iP(At) = e" iAt9 ^^lV(0). (2.4) 
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This is because (|2.3[) exactly conserves \ip\ 2 (multiply (]2.3[) by ip*, the complex conjugated equation by if> and subtract) 
and the nonlinear potential is just a constant in (|2.3p . This is also clear from 



hHA£)| 2 = |VK0)| 2 , (2.5) 

since ip(0) is only multiplied by a phase. Eq. (|2.2j) and (|2.4[) are the basic building blocks for constructing splitting 
algorithms for solving the nonlinear Schrodinger equation. Eq. (|2.4p is the fundamental justification for using the 
"latest" wave function in computing the nonlinear potential 3 . (See also below). Define a time-independent operator 
V such that 

v\m) = 9\m\ 2 \m)- m 

Note that V only acts on \ip(t)} and does not act on its own eigenvalue g\ip(t)\ 2 . It follows that 

e- lAt9 \iP(t)) =e- iAt9 ^W 2 \i>(t)). (2.7) 

The crucial point here is that V has no time-dependence, when it acts on any \ip(t)), it produce the eigenvalue 
g\ip(t)\ 2 . The resulting time-dependence of the nonlinear potential is due entirely to the state vector \ip(t)) and not 
to the operator V. The exact solution can then be written in operator form as 

|^))= e -^ f +^(0)>. (2.8) 

For our purpose here, we only need to know (|2.7p and not the explicit form of V. For an elegant, but rather abstract 
construction of V, see Strauch's^ recent work. 

III. DERIVING SPLITTING ALGORITHMS 

To solve (|2.8|) by splitting algorithms, one factorizes the evolution operator to any order with a suitable set of 
coefficients {U,Vi} via 

e e(f+V) = TJ e ^ eE ^ (31) 

t 

where we have denoted e = —iAt. For example, we can have the second order algorithm 2A as 

= e^e ef e^V(O) = e^l^V f >(0) (3.2) 

where according to ()2.4j) or (|2.7p . we must take 

= c^e^^ 2 iP(O). (3.3) 

Algorithm 2A only requires one-pair of FFT (forward and backward) to achieve second-order accuracy, which is the 
same number of FFT needed for a first-order algorithm. If the nonlinear potential is treated as a time-dependent 
potential, as done by Bandrauk and Shen^, then we would have the algorithm^ 

^(At) = e -^A t V(Ai) e -iAtf e -i|A t y(O)^ (0)> ( g 4) 

In this case, since the last factor is only a phase, 

V(At)=g\ib(At)\ 2 =g\4>\ 2 , (3.5) 

the result is the same as (|3.2[) . If one ignores the time-dependence^ and uses V(At) = V(0) = .g|-0(O)| 2 , then algorithm 
(|3.4p is degraded to first order. 

Similarly one has the second-order algorithm 2B, 

i/>(£t) = e^ et e eir ^ ef Vj(0) = e^ ef e £9l<#,|2 e^ ef V(0), (3.6) 

where here 

$ = e i ef ip(0). (3.7) 
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In the time-dependent potential approach, one would have instead, 

if>(At) = e~^ Atf e -^tv(At/2) e -<^Atf ^ (Q) _ (3 8) 

One must now evaluate V(At/2) = g|-i/>(Ai/2)| 2 . Since the algorithm is only second order, one can simply approximate 
the midpoint wave function to first order, 

^(At/2) = e-' l i Atv( - At ^e~ l i Atf ^(0), (3.9) 

and therefore 

|V>(At/2)| 2 = |e-^ Atf i/;(0)| 2 . (3.10) 

Again, the result is the same as (|3.6j) 

For fourth and higher order algorithms, the time-dependent potential approach cannot be easily implemented. It 
is much more efficient to use the "latest" intermediate wave function than to estimate the intermediate-time wave 
function to third or higher order. Thus higher order algorithms are currently possible only with the use of the 
time-independent formalism based on the original finding of Taha and Ablowitz. 

The fourth-order Forest-Ruth (FR) 9 algorithm, which is the triplet concatenation ^ 10 ' 11 of algorithm 2A 

T FR {e) =T 2A { Cl e)T 2A {cQe)T 2A {cie) (3.11) 

with ci = 1/(2 — 2 1 / 3 ) and cq = — 2 1 / 3 /(2 — 2 1 / 3 ) has been verified by Javanainen and Ruostekoski as obeying the 
"latest" intermediate wave function rule. However, since this triplet concatenation will convert any second-order split 
algorithm to fourth-order, verifying this algorithm alone does not constitute a check on more general fourth-order 
algorithms. (Recall that algorithm 2A can also be derived from the time-dependent approach without explicitly 
invoking the "latest" wave function rule.) (Javanainen and Ruostekoski have also verified the "latest" wave function 
rule on a class of third-order algorithms independent of 2A.) To seal this loop-hole in our verification process, we also 
consider more general fourth-order algorithms previously studied by McLachlan— with 9 operators, 

T M = ... exp(et V r ) exp(ewif) exp(rtiy) exp(ev 2 f) exp(et 2 V). (3.12) 

The factorization is left-right symmetric and only operators from the center to the right are indicated. The fourth-order 
order condition requires^ 3 - that 

vi = ^-v 2 , i 2 = i-4Mi, t = l-2(t 1 +t 2 ), (3.13) 



1 / /9ti -4±2w 



W = ^j3-m 1+ 9tl V2 = -{1 T] J^— j (3.14) 

and that the free parameter t\ < 0. This algorithm requires 4 pairs of FFT but has a much smaller energy error and 
greater stability than that of FR. (The coefficient designation does not match the the operators because the algorithm 
has been adapted from its classical version by interchanging T <-* V .) There are four solution branches for v 2 . The 
choice of 

121 , 

<1 = 3924 (12 - V ^ ) "-°- 299 

with 



1 / 9t 1 -4 + 2w . 
"2 = 7 ! + \/ 51 ( 3 - 15 ) 



4 V V 3ii 

reproduces McLachlan'si^ recommended algorithm. By varying t\ and using different branches of v 2 , it is possible to 
optimize the algorithm for specific applications. For application in the next section, the results are not very sensitive 
to the branch of v 2 nor the choice of t\, as long as t\ is in the range of [-0.1, -0.4]. More higher-order splitting 
algorithms can be found in Refsi 14 ' 15 ' 16 ' 17 . 
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IV. NUMERICAL VERIFICATIONS 

To verify the order of convergence of these algorithms, we apply them to the Gross-Pitaevskii equation with a 
harmonic trap in ID, 

To gauge the accuracy of any algorithm, we monitor the fluctuation of the total E, 

r 00 1 rl 2 1 1 

E = J J x r(t)(--^ + -u J ^ + -gm)\ 2 )m (4.2) 

If the time evolution is exact, E would remain a constant. For u = 1, g = 5 and ip(0) = tpo(x), the ground state wave 
function of the harmonic trap, the initital total energy is 

E=- + — = « 1.497355701. (4.3) 
2 2V2^ 

The x-interval used is [-20:20] with 2 9 = 512 grid-points. The results are unchanged if one doubles the grid-points. In 
Fig.l we plot E as a function of time for algorithm 2A at At — 0.05 and At = 0.025 . One observes that the energy 
fluctuation at At = 0.025 is about 1/4 of that at At = 0.05, as befitting a second order algorithm. The results for 
fourth-order algorithms FR (Forest-Ruth) and M (McLachlan) at At = 0.05 are also shown. It is clear that even if 
one take 1/4 of algorithm 2A's error at At = 0.025, corresponding to At = 0.0125, that error is still much larger than 
those of fourth-order algorithm FR and M (i.e., running algorithm 2A four times at At = 0.0125, using 4 pairs of 
FFT, would still be inferior to algorithm FR which uses only 3 pairs of FFT). 

In Fig. 2 we greatly magnified the scale so that the fluctuations in the fourth-order algorithms are also visible. This 
time, when the step size of algorithm FR is half, the error in E is reduced by a factor of 16, confirming the fourth-order 
convergence of the algorithm. The energy error of algorithm M at At — 0.025 is ~ 10 -6 , which is too small for a 
visual comparison. 

In both Figures 1 and 2, the total energy eventually blows up for all calculations, despite the fact that total energy 
error is only 10~ 6 for McLachlan's algorithm. This instability is directly related to the strength of the nonlinear 
potential. The rather large value of g = 5 was chosen so that the instability would show up after a short run. The 
energy blow up can be delayed, but not eliminated, by reducing At. (See further discussion in Section VI.) 



V. THE CAUSE OF INSTABILITY 



The eventual instability as shown in Figures 1 and 2 demands an understanding of its fundamental cause. To 
study this, we decompose the general wave function into Fourier components and focus on the propagation of a single 
component with wave vector p in ID, 



i>(x,t)=Ae ipx - iu,t . (5.1) 



This is a solution to (|2.ip if u> is given by 



uj = -p 2 +g\A\ 2 =E p + U, (5.2) 



where we have denoted E p = ^p 2 and U = g\A\ 2 . Suppose now the spatial part of ip is contaiminated, due to 
numerical errors, by very small amplitude, side-band wave vectors p + k and p — k so that 

iP{x) = Ae ipx + ae i( - p+k '> x + 6e l(p ~ fe)a; , (5.3) 

how will the error amplitudes a and b be propagated by splitting algorithms? (This side-band analysis was inspired 
by the classical work on Fourier analysis of nonlinearly interacting waves^.) The effect of e _zAtT on ip(x) is trivial; 
all amplitudes are multiplied by a phase, 

A' = e- lAtE ?A 

a' = e~ lAtE "+ k a (5.4) 
V = e- lAtE ?- k b. 



To compute e lAtv ip(x), one must compute |^>(x)| 2 using (|5.3p . The result, by keeping terms only to first order in a 
and 6, is 

A' = e~ tAtu A 

a' = e~ tAtu [a-iAt(Ua + gA 2 b*)} (5.5) 
b' = c- tAtu [b-iAt{Ub + gA 2 a*)]. 

Thus the first order splitting algorithm e~ lAtv e~ lAtT ip(x) modifies the amplitudes by composing (15. 4|) with (|5.5p . 
yielding 

A n+1 = e-^ E ^A n (5.6) 

a n+1 = e -i*KE p+k -E k +u)[ ane -iAtE k -iAtU(a n e- lAtEk + {b n e' lAtEk )*e' as -)] (5.7) 

b n+1 = e -i±t{E p - k -E k +u)[ bne -iAtE k -iAtU(b n e- lAtEk + {a n e~ lAtEk )*e~ l25 ")] 1 (5.8) 

where we have defined 

A n = \A n \e~ t5 ". (5.9) 
The algorithm correctly propagates A and preserves the norm |j4|, 

A n = e -inAt(E p +U) Ao _ (5.10) 

For notational clarity, we will take Aq to be real with = so that we don't have to keep track of this initial phase, 
yielding 

S n = nAt(E p + U). (5.11) 

(Keeping the initial phase simply transfers it to subsequent amplitudes and has no bearing on the issue of instability.) 
To see the growth in a and b, we factor out their overall phases as follow 

_ ~inAt(E p+k -E k +U) 

b n = e -inAt(E p _ k -E k +U) ( 5 12 ) 

and reduce l|5.7j) and (|5.8[) to 

a n+1 = a n e- lAtEk -iAW{a n e- lAtEk + {(3 n e- lAtEk )*) (5.13) 

Pn+i = P n e- lAtE » ~ iAW(f3 n e~ lAtE * + (a n e^ AtE ")*). (5.14) 



These two equations can also be interpreted as a first-order splitting algorithm, with the "kinetic" term giving 

a' = e~ tAtE *a 

P = e~ lAtEk P (5.15) 



and the "potential" term producing 

a = a~iAtU(a + 8*) 

3' = 8-iAtU(f3 + a*). (5.16) 
A closer examination reveals that (|5.15|) and (|5. 16|) are exact solutions to following equations 

da dB 

i^=E k a, i-£ = E k 0, (5.17) 

da dB 
i^L =U(a + P*), i-- = U(8 + a*). (5.18) 
dt dt 

Thus the algorithm is trying to solve the original unsplitted equations 

df\ 

i^ = (E k + U)a + UB* 
dB 

i-£ = (E k + U)B + Ua*, (5.19) 
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which have general solutions of the form 

a = cc- lUkt + de lUkt , (5.20) 

with 



Q k = y/E k (E k + 2U). (5.21) 

This is the famous Bogoliubov spectrum^ of elementary excitations in a uniform Bose gas. It shows up here be- 
cause the nonlinear Schrodinger equation is just the Gross-Pitaevskii equation for describing a uniform Bose-Einstein 
condensate 2 ^. The Bogoliubov spectrum in the current context, is the background "noise" excitations of the nonlinear 
Schrodinger equation. If one were able to solve (|5.19[) exactly via (|5.20p . there would be no instability because the 
amplitude of a in (|5.20p is finite. However, when (|5.19p is solved by splitting, (|5.16p no longer preserves the norm 
and the modulus of these error terms at selected ranges of k will grow exponentially. 
To study this growth, take (3q = ao, so that the splitting forms (|5.15|) (|5.16|) simplify to 

a' = e' lAtEk a (5.22) 
a' = a-iAtU(a + a*). (5.23) 

Now we assert without giving a detail proof that beyond first-order, for any splitting algorithm in solving the nonlin- 
ear Schrodinger equation, the error Fourier components will grow correspondingly according to splitting (|5. 22115. 23)) 
with the same splitting coefficients. For example, corresponding to algorithm 2A, the growth of the error Fourier 
components is given by 

ai = a - i^AtU(a + a* Q ) 

a 2 = e - lAtBfe ai (5.24) 
a 3 = a 2 - i-AtU(ct2 + a* 2 ) 

The subscripts here simply label the individual steps in the algorithm. The last labelled value is the updated variable 
after one time step. Denoting this updating as £2A(At), the error growth of the Forest-Ruth algorithm is then 

£ FR {At) = E 2A {c 1 At)E 2A {c At)E 2A {c 1 At) (5.25) 

and McLachlan's algorithm as 

a\ = a - i(t 2 At)U(a + ckq) 

a 2 = e-W^ax 

«3 = ct2 — i{t\At)U [a 2 + a 2 ) 

—iv-tAtE k 

, etc. (5.26) 



To verify the validity of our assertion, we run the normal algorithm on an initial wave function having the p = 
component with amplitude A = 1, and all other Fourier components set to e -25 , at g = 5 and At = 0.2. The 
resulting Fourier amplitudes are then outputted every time steps for seven time steps. Their modulus are shown as 
plus signs for the above three algorithms in Figs. 3-5. Instead of plotting the magnitude of these Fourier amplitudes 
as a function of fc, we plot them as a function of AtE^/it, which is more revealing. Also plotted as solid lines, are the 
predicted error amplitudes given by ()5.24|) . (|5.25j) and (I5.26|) for seven time steps. The perfect agreement in all three 
cases confirms our assertion and our side-band analysis. 

To understand the pattern of instability as shown in Figs. 3-5, we rewrite the splitting forms (|5.22p and (I5.23P as 
matrices acting on the real and imaginary part of a 



T K a ^)< («?)= v(a<) (")' (5 - 27) 



with 



T(At) = f _ C s S c ) , V(At) = ( _\ u J ) , (5.28) 
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and where we have defined 



c = cos(x), s = sin(ai), x — AtE k , and u = AtU. 
The updating matrix corresponding to algorithm 2A is therefore 

M 2A (At) = V(±At)T(At)V(±At) 

c — us s 
(u 2 — l)s — 2uc c — us 



(5.29) 



(5.30) 



This is a special form of a matrix with equal diagonal elements and unit determinant. This is due to the left-right 
symmetric form of the matrix product ( i.e., the algorithm is time-reversible2i) and that both T and V have unit 
determinant. Such a matrix has the special property that its eigenvalue is given by 



ei, 2 = C ± VC^l, 



(5.31) 



where C is just the diagonal element (or half of the trace of the matrix). If |C| < 1, the eigenvalues are complex 
with unit modulus and the algorithm is stable. If |C| > 1, the eigenvalues are real with one eigenvalue always greater 
than unity. Thus by just plotting C against x — AtEk, one can immediately determine the regions of instability. For 
algorithm 2A, we have 



C(x) = cos(x) — u sin(a;) = Co cos(x + 8). 



with 



Co = yl + u 2 and 8 = tan 



(5.32) 



(5.33) 



It is then immediately clear that as long as u ^ 0, the algorithm is unstable for x in the interval [nn — 28, nir] where 
n = 1,2,3... etc.. At a fixed U , decreasing At reduces u and 8, and hence the width of the instability region, but does 
not remove the instability (but see further discussion in the next section). In Fig. 3, this C-function is plotted and 
lowered to -28 so that the interval where |C(x)| > 1 can be directly compared with the observed regions of instability. 
The peak instability occurs at x = nix — 8 with the maximum eigenvalue 



d,2 



vT 



(5.34) 



For At = 0.2 and U = 5, we have u = 1, 5 = 7r/4 and |e| = 1 + V%- After seven iterations, the e-fold increase of the 
peaks would be log((l + v2) 7 ) = 6.16962, which is the six e-fold increase of amplitude observed in Fig. 3. Thus we 
have completely accounted for, both qualitatively and quantitatively, the pattern of instability as shown in Fig. 3. The 
corresponding C-functions for the Forest-Ruth and the McLachlan algorithm are also plotted in Fig. 4 and 5. Their 
C-functions are too lengthy for a written display. (The analytical expression for McLachlan's C-function is more than 
a page long using Mathematica.) 

By comparing Fig. 3 and 4, one sees that the Forest-Ruth algorithm has a greater error growing rate than 2A. We 
will see in the next section that this is precisely the reason why the FR algorithm blew up earlier than 2A in Fig.l. 
Finally, as shown in Fig. 5, McLachlan's algorithm manages to shift the C-function is such a way that the error peaks 
at x/ir = 1,3 are nearly eliminated. 

Further insights into the origin of this instability can be gained by representing T(At) and V(At) in terms of 
traceless matrices, 



T(Ai) = exp 



At 



E k 
-E k 



V(At) = exp 



At 





-2E7 



One can then immediately identify the unsplitted evolution operator as 



exp 



At 







-E k - 2U 



cos(f2 fe Ai) sm(n k At) 
-sm(n k At) cos(f2 fe Ai) 



(5.35) 



(5.36) 



which is that of a harmonic oscillator with the Bogoliubov spectrum fl k - Were one able to split it alternatively as 



T'(Ai) = exp 



At 



E k 




V'(At) = exp 



At 











-E k -2U 



(5.37) 
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one would recover the stability criterion normally associated with the harmonic oscillator. For example, the corre- 
sponding second-order algorithm 2A, V'(iAi)T'(Ai)V'( i Ai), would then yield a C-function of 

C = l-ifi|Ai 2 , (5.38) 

which limits stability to At < 2/fifc, a well known result. This limit is actually worse than x < tt — 28, which, as 
U — > 0, is At < ir/E k . Our original splitting (|5.35p is therefore better the usual harmonic oscillator splitting (|5.37[) . 
Moreover, in contrast to Fig. 3, the usual harmonic oscillator splitting would have no stable region whatsoever beyond 
AtE k >n\ 

In this section we have shown that the error growing pattern of any splitting algorithms when solving the nonlinear 
Schrodinger can be analytically understood. The instability is due to the exponential amplification of high k noises 
at E k >n/At. 

VI. THE INSTABILITY OF THE GROSS-PITAEVSKII WAVE FUNCTION 

We now repeat the calculations of Fig.l at At = 0.05 for 1200 time steps, to the point where the algorithm FR 
begins to blow up. We plot in Figs. 6-8, the modulus of the fc-space wave function 1^(^)1 as a function of AtE k /ir at 
every 100th time step. The initial Gaussian wave function is the straightline seen plunging down close to vertical axis. 
Because of limited numerical precision, that line levels off to some random values around e -35 » 10~ 16 at high E k . 
These are the initial random errors of the wave function. When the algorithm acts on the wave funtion, these random 
errors are amplified successively and grow in time. For algorithm 2A, Fig. 6 shows error peaks at x/tt = 1,2 and 4, 
which is in agreemwnt with Fig. 3, but no discernable peak is seen near x/tt = 3. For the Forest-Ruth algorithm, 
Fig. 7 shows a promenient peak at x/tt = 1, followed by a peak-shoulder structure at x/tt = 2 and 4, in agreement 
with Fig. 4. For McLachlan's agorithm, Fig. 8 shows that the error peak at x/tt = 1 is conspicuously absent, and only 
peaks at x/tt — 2,4 are visible. This is in excellent agreement with the predicted error structure of Fig. 5. In the case 
of the Forest-Ruth algorithm, the error peak at x/tt = 1 has grown sufficiently to distort the wave function and cause 
the energy to blow up. These exponentially growing error peaks are like ticking time bombs, harmless at first, but 
eventually overwhelm and destroy the wave function. 

For a continuum wave function, this instability is unavoidable as long as At is finite. However, for a discrete wave 
function defined at only N grid points, there is a loop-hole. For a finite A-point calculation, the maximum k vector 
is k max = Ntt/L so that AtE k /iT extends only out to (0.05)0.5(5127r/40) 2 /7r ~ 12.9, as shown in Figs. 6-8. Thus one 
can take advantage of this and force stability by making At so small that 

AtE™ ax < x mm , (6.1) 

where x m i n is the smallest value of x such that |C(x)| = 1 and E™ ax = \k r ^ nax . For most algorithms at small At, 
Xmin ~ tt- This criterion (|6.ip simply shrinks the entire range of E k values to below the first instability point. Thus 
the RF calculation would be stable for At < 7r/(0.5(5127r/40) 2 ) = 0.0039. A more refined calculation at higher N 
would required an even smaller At. Such as small At would make long-time simulation very time consuming. On 
the other hand, (|6.ip also implies that stability can be achieved by lowering k max , i.e., using fewer grid points. For 
example, at N = 128, 7r/(0.5(1287r/40) 2 ) = 0.062. When the FR algorithm is rerun at At = 0.05 but with N = 128, 
the total energy is indeed stable out to t = 300. However, the wave function now looked very jagged. Thus for long 
time simulations, one muct choose At and N judiciously. 

The instability observed here is very similar to the "resonance" instability of multiple-time step algorithms used 
in biomolecular simulations 22 . There, stability requires that At < tt/uj, where u> is the faster physical frequency in 
the problem. The latency in the energy blow-up has also been observed in density functional calculations using split 
algorithms 23 . The energy blow-up there is more gradual, but it is undoubtedly related to the nonlinear Kohn-Sham 
density used, for which the nonlinear Schrodinger equation is the simplest prototype. 

VII. CONCLUSIONS 

In this work we have shown how splitting algorithms of any order can be devised to solve the nonlinear Schrodinger 
equation. The key ingredient is the exact solution of the potential equation (|2.4[) . as pointed out earlier by Taha 
and Ablowitz 1 . This explains Javanainen and Ruostekoski's finding 3 without the need to construct Strauch's special 
operator 4 . Solution (|2.4p clearly generalize to the case where g\ip\ 2 — * v(\ip\), implying that this class of general 
nonlinear equations can also be solved by splitting algorithms. 
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In the course of verifying these alogrithms by solving the Gross-Pitaevskii equation, a latent instability is observed 
in all the algorithms. This instability persists regardless of the order of the algorithm and despite excellent total 
energy conservation. A detail error analysis reveals that this instability is intrinsic to splitting algorithms and can 
only be avoided if (|6.1[) is satisfied. 

The main advantage of higher-order algorithms is that a larger At can be used for more efficient simulations. 
However the stability criterion (|6.1|) dictates a small At regardless of order, thus negating much of the presumed 
advantage of using higher order algorithms. (Of course, higher order algorithm are useful for short time simulations, 
where results can be obtained prior to the blow-up.) This work also suggests that one must not use just any higher 
order algorithm, such FR, but higher order algorithm with a higher x m i n , such as McLachlan's algorithm. How 
algorithms can be derived systematically with a higher x m i n is a fitting subject for a future study. 
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FIG. 1: (Color online) The fluctuation in the total energy E (|4.2|) . when solving the real time Gross-Pitaevskii equation by 
second order algorithm 2A and fourth-order algorithms FR (Forest-Ruth) and M (McLachhan). 
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FIG. 2: (Color online) A magnified view of the fluctuation in the total energy of two fourth-order algorithms FR and M at two 
time-step sizes. 
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FIG. 3: (Color online) The growth of the error Fourier amplitudes due to algorithm 2A for seven time steps at g — 5 and 
At = 0.2. The plus signs denotes the algorithm's actual output; the seven solid lines are the predicted error from the side- 
band analysis (|5.24|l for seven time steps. Centered on -28 is the algorithm's C-function for predicting regions of stability and 
instability. 
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FIG. 4: (Color online) Same as Fig. 3 but for the Forest-Ruth algorithm. 



The predicted error is given by (|5,25|) . 




FIG. 5: (Color online) Same as Fig. 3 but for McLachlan's algorithm. The predicted error is given by (|5.26|) . 
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FIG. 6: (Color online) The modulus of the Gross-Pitaevskh momentum wave function |^>(fc)| at every 100th time-step due to 
algorithm 2A. The time step size is At = 0.05. 



17 




_30 I l i i i i i i ^ I 

2 4 6 8 10 12 14 

AtE k /7l 



FIG. 7: (Color online) The modulus of the Gross-Pitaevskii momentum wave function |^>(fc)| at every 100th time-step due to 
Forest-Ruth algorithm. This is the momentum wave function corresponding to the energy calculation of Fig.l up to t = 60. 
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FIG. 8: (Color online) Same as Fig. 7 for McLachlan's algorithm. 



